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SUMMARY 

An analysis for modeling damolng seals for use In Space Shuttle main engine 
turbornachlnery Is being produced. Development of a computational fluid mechanics 
code for turbulent, Incompressible flow is required. This document Is an interim report 
on the first year of work. It wos prepared at Continuum, Inc. for NASA-MSFC under 
Contract NAS8-35508. 
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INTRODUCTION 

One of the limiting factors In the operation of the Space Shuttle main engine (SSME) Is 
the speed of the turbopumps The turbopumps currently operate near the whirl Instabllty 
speed. Further Increases In S5ME thrust will require an increase In this speed. It has 
recently been suggested that the creation of turbulence In the leakage flow of the 
current seals would provide additional damping for the rotor. This Idea was Investigated 
analytically in Ref. I. The current work is to use numerical methods to model these 
damping seals and thus provide a tool to further analyze this concept. 

In order to describe the flow in a damping seal, a variational method (known as the VA5T 
code, VAriational Solution of Transport equations) will be used to solve the governing 
equations. The VAST code was developed at Continuum, Inc. based on the work of 
Prozan (Refs. 2,3). For this work D rozan solved the compressible form of the equations, 
and used the Ideal gas equation of state. The investigation of the damping seals requires 
an incompressible flow solver with liquid as the working fluid. 

Incompressible flow solution codes are almost universally implicit in nature. In most 
codes pressure is generally obtained by solving one of the momentum equations, or by a 
Poisson equation. Attompls to develop explicit pressure solutions are primarily 
characterized by the artificial compressibility of Chorin (Ref. 4). Recent applications of 
this method to the SSME have been demonstrated (Refs. 5,6). 

In developing an incompressible version of the VAST code, a pressure solution which falls 
within the framework of the variational statement is required. A first attempt to 
include an explicit pressure solution in VAST was made by Farmer (Ref. 7). This 
procedure utilized one of the momentum equations to basically treat gravity flows, but 
this causes the code to be directionally biased. This method generally gives different 
answers if different paths of integration are used (Ref. 8). 

In describing the damping seal, the pressure is one of the most critical variables, whereas 
many of the aforementioned codes use pressure simply as a means of arriving at 
velocities. Indeed, some researchers have eliminated pressure entirely in their 
incompressible solutions, as in Ref. 9 or as with the stream function - vorticity 
formulation. 
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In the following work, efforts at obtaining a pressure solution that fits the requirements 
of the damping seal model, and the subsequent development and evaluation of the VAST 
code for Incompressible flows Is described. Since the flow In the seals Is turbulent, 
turbulence models must also be Included In the code. Thus, work supporting turbulent 
boundary layer modeling Is also discussed. Some results are presented and remaining 
work Is outlined. 


THEORY 

In the current approach to obtain a solution to the Incompressible problem, the general 
form of the compressible Navier-Stokes equations are initially considered. Thus, the 
continuity and momentum equations are, respectively: 



for a single species and neglecting body forces. The energy equation can be used if 
temperature changes are desired. However, since an incompressible code is sought, the 
energy equation decouples from Eqs. (I) and (2) so that this extension is 
straightforward. For simplicity, only the isothermal case is considered here. 

In the usual form for incompressible flow, Eq. (I) Is written as 



( 3 ) 
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This equation implies an Infinite wave propagation speed. Consequently, any type of 
truly unsteady, explicit solution to the governing system of equations Is precluded. 
Because of the nature of the damping seals problem, It Is desirable to obtain a set of 
equations that will allow unsteady analyses. 

Pressure Equations 


Several methods for obtaining accurate pressure for incompressible problems have been 
studied. One of the most widely used methods Is the Poisson equation (Ref. 8). The 
differential form of this equation Is 


V = 


(4) 


where Sp |$ a function of the velocities as well as the first three derivatives of the 
velocities. Reference 8 provides a detailed presentation of Eq. (4), along with stablity 
considerations and boundary conditions. Eq. (4) can only be solved Implicitly, although It 
may be combined with explicit solutions of Eq. (2) and used iteratively to obtain some 
time varying solutions. 

Another solution procedure which has been examined is the penalty function method. An 
excellent review of this method is given in Ref. 10; however, this procedure is also 
implicit, and in its usual form is a true finite element formulation. The variational 
procedure, on the other hand, consists of both finite element and finite difference 
characteristics, with the actual time integration being carried out with differencing 
schemes. 


A new approach which was considered for this problem was the use of mechanical 
energy. Since minimum kinetic energy is a known variational principle for some special 
cases of incompressible flow (Refs. 9,11), it was hoped that this approach could be 
utilized. The mechanical energy is: 


IJ Iff 


3 ( . 5pV z ) 


du - 


3 t 


JJ\(- 5pV 2 )V - x • V j * d s (5) 
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As with the previously discussed methods, It Is seen that there Is no time-dependent 
pressure term once again dictating an Implicit solution. Furthermore, Eq. (5) does not 
actually provide new Information beyond Eq. (2), and there is no a priori satisfaction of 
continuity, while Eq. (4) does satisfy continuity In the form of Eq. (3). 

One of the earliest attempts to develop explicit pressure equations for incompressible 
codes was the artificial compressibility of Chorin (Ref. 4). This work was restricted to 
steady problems, and had the drawback that the assumed compressibility had no basis In 
physical reasoning, but rather played the role of a relaxation parameter. Recent work by 
Chang and Kwak (Ref. 5) and Kwak, et al (Ref. 6 ) have furthered the applicabllty of this 
procedure. The continuity equation can be obtained in the form 



du + 8 


IP - 


0 


( 6 ) 


The explicit nature of the equation is readily apparent, and 8 is the artificial 
compressiblity coefficient. Again, 8 is a term introduced for convenience, although in 
Ref. 5 some theoretical bounds are developed by using an artificial speed of sound and by 
introducing artificial waves. These bounds are also dependent upon the type of numerical 
differencing scheme used (Ref. 5). Nonetheless, the procedure was applied successfully 
to some significant flow problems, including some 5SME applications. An attempt to 
develop a time-accurate solution was also made in Ref. 6 by replacing zero on the right 
hand side of Eq. (6) with another pressure term and Iterating. This procedure, however, 
is still in the developmental stage. 

Although the form and ease of application of Eq. (6) is very appealing, the lack of 
physical rec ->ning to obtain the equation, as well as the trial and error procedure often 
necessary to determine an acceptable value of 8 , detracts from Its usefulness. This 
deficiency has been overcome by the work of Kawahara and Hirano (Ref. 12). This 
analysis will not be presented here, but rather a very similar apprach which parallels 
their work is presented. Differences in the resulting pressure equation will then be 
pointed out. 


- 4 - 
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The initial step In the derivation Is to obtain an equation of state. The Idea of a limited 
compresslbllty Is used and the assumed functional form of the equation of stare Is 


P = Pip) 


(7) 


For a liquid (or gas, also), the bulk modulus of elasticity expresses the compresslbllty of 
the fluid, 


dP 

E = 

P 

The finite acoustic speed In a liquid is given by 


P 


(8) 


(9) 


Hence, for all fluids with a finite acoustic speed (which In actuality encompasses all 
fluids), Eqs. (7), (8), and (9) may be used to write the equation of state as 


dP = c 2 d p 


( 10 ) 


For most of the following work, Eq. (10) is used in the differential form wherever 
possible. However, In some instances it Is necessary to integrate this equation with 
respect to some reference state. This reference state, P r and p , must be taken to be 
the usual theoretical Incompressible values of P and p for the fluid being considered. 

If the Integrated form of Eq. (10) is substituted into the continuity equation, Eq. (I), the 
result, after some manipulation, is 


- 5 - 
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Iff I 7 du * // F " V - 


ds + ( p 


- P 


■’ IP ■ 


(is S 0 


(II) 


Note that In this equation, the only apearance of density is the reference density p r . 
Eq. (I I) Is a physically meaningful explicit pressure equation. If the quantity 


(P-P r ) 




( 12 ) 


Is subtracted from the left hand side of Eq, (II), the resulting equation will be Identical 

to the method of Ref. 12. Equation (II) Is also very similar to the artificial 

2 

compressibility equation. If 8 of Eq. (6) is replaced by P r c - P r and the 

PV • ds term Is neglected, the equations become the same. As with artificial 
compressibility, It is also noted that as the acoustic speed, c, approaches infinity, Eq. 
(II) reduces to the incompressible continuity equation, Eq (3). 

In solving for the momentum, Eqs. (2), the incompressible density, p , Is directly 
substituted for p . A first order analysis indicates that this assumption is valid, 
however, more study is required. 


Turbulence Modeling 

The primary direction of work for turbulent flow modeling has been the development of 
shear stress wall functions for use in the VAST code which will permit viscous analyses 
without requiring an excessive number of computational grid points. It is impractical to 
resolve the flowfield in the vicinity of the wall with enough grid points to accurately 
calculate the wall friction. Wail functions are commonly used to provide the required 
boundary conditions; however, care must be exercised In order not to make the wall 
functions too empirical. The end result is to predict frictional losses from the detailed 
flow vectors, not from the mean channel flow. The following procedure was developed 
and tested wi+h geometrically simple problems to provide the necessary CFD tools for 
SSME analysis. 
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For turbulent wall flows for a smooth wall, Including fully developed pipe flows, the 
following empirical velocity profiles are vatldr 


u + = 5.5 + 2.5 !n(y + ) for y + > 30 


(13) 


u + = -3.05 + 5.0 I n ( y + ) for 5 < y + < 30 (14) 


u + = y + for y + < 5 

where 

+ / * 

u = U ' U 

. tx 0 /p )°- 5 

y* = (t 0 /o)°- S y/v 


( 15 ) 


Blasius* empirical shear stress relationship is appropriate. 

t q = 0.0255 pu 2 (v/u^)°‘ 25 (16) 

where l - R for pipe flow 

1= 6 for boundary layers 

The boundary layer thickness Implied by Eq. (16) and a 1/7 power-law profile is 

6 = 0.376 X/R^° * 2 (17) 

From Eq. (15) 


9 u 


( — ) 

3 y 
'w 


T 0 /p = (uV/v 


(18) 
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In terms of real distance from the wall, Eq. (13) represents most of the boundary layer, 
therefore the following computational procedure Is sugg<?3led. A flctltous wall Is 
assumed to be 0.0005 ft away from the real wall, and no flow Is assumed to occur 
between the two walls. Equation (13) Is valid at .0005 ft from the waltj hence, If u Is 
calculated with a slip boundary condition, u* Is determined. Equation (18) Is used to 
calculate the velocity gradient at the wall. Since Eq. (13) Is not explicit In u , the 
approximation 


0. 1662529 O 0*867325 /(y/v) 0. 132675 


(19) 


Is used. These equations determine the velocity gradient and shear stress at the wall. 
An eddy viscosity Is used to determine both the local sheer stress and the variation of 
this stress with distance from vhe wall, y. 


P T = 0.07 u *£ p ( FR ) + p (20) 

where 

FR = (y/0.3 £ ) for 0 « (y /£) <,0.3 
FR = 1.0 for ( y / £ ) >0.3 

Equations (19) and (20) and the momentum equations were used to calculate the turbu'ent 
boundary layer over a flat plate between one and two ft running length over the plate. A 
turbulent boundary layer was assumed at the leading edge of the plate. The flow was air 
with a free steam velocity of 100 fps (this Is an approximate Reynolds Number of 10^). 
By adjusting the constants in Eqs. (19) and (20), the profile at the end of the plate was 
predicted to be that shown in Fig. I. The fit is very good, especially near the wall; the 
calculated wall shear stress Is within 5% of the correct answer. This procedure Is 
accurate enough to extend its development to more geometrically complex flows. The 
reasons for the necessity of adjusting the constants in Eqs. (19) and (20) and for the lack 
of better fit at y's near the free stream side of the boundary layer are still under 
Investigation. 


- 8 - 


Velocity (fps) 



V (Ft.) 


Fig. 1 Turbulent Boundary Layer Calculation 
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Solution Procedure 


In order to use the variational approach developed by Prozan (Refs. 2,3), the governing 
equations, Including the pressure equation must be obtainable from g fractional. 
Functionals are already known for some special cases of Incompressible flow (Ref. 1 1), 
generally employing the kinetic energy. However, since limited compressibility Is used In 
at least two of the pressure equations previously developed, the entropy functional was 
used to retain as much commonality as possible with the earlier developments by Prozan. 

Beginning with the entropy equation 


dS = - dE + PdV 
T 


( 21 ) 


the differential form of the equation of state, Eq. (10), may be used and, after 
Integration, the entropy functional Is obtained as 


PS 


c I n T 
v 




„ 2 , , 2 
pc I n ( c - 


P 

_) 

P 


+ const. 


( 22 ) 


Now, P, u and v are the independent variables, and partial derivatives of Eq. (22) with 
respect to each of the independent variables are needed. For the velocity derivatives, a 
velocity-pressure relation Is required, and Bernoulli's incompressible relation Is assumed. 


The derivatives are then derived as 


3(pS) 
3 P 

3 (pS) 
3 (pu) 

3 (pS ) 
3 (pv) 


c y I n T I n ( c - _) 


u 

T 


I - 


c 2 + . 5 V 2 


v r 

T 


c 2 + . 5 V 2 


(23) 

(23) 
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Although Eq, (22) is developed specifically for the pressure equation as given by Eq. (I t), 
the analogous development for the artificial compressibility scheme, Eq. (6), is mado by 
simply using th* equation of state developed for this method. (See Refs, 4 and 5 for this 
equation of state). 

The system of equations, Eq. (2), coupled with the preferred pressure equation, can now 
be solved numerically with the use of proper boundary conditions. Boundary conditions 
will be discussed In more detail In the next section. Equations (23) are used to prescribe 
the numercial scheme when Eq. (I I) is used for the pressure. Similar procedures are used 
when alternate pressure equations are employed. 


MODEL EVALUATION 

The various pressure solutions discussed in the previous section were evaluated for 
possible use in the damping seals analysis. The mechanical energy equation, £q. (5) and 
penalty function methods were eliminated early, primarily due to the implicit nature of 
the equations. Mechanical energy also suffers from other deficiencies as previously 
pointed out, particularly with regard to the satisfaction of continuity. The most 
promising approaches are the Poisson equation and the use of artificial or limited 
compressibility. 

The use of a Poisson equation to obtain pressures is widely known and numerous 
applications with different numerical codes are thoroughly discussed in Ref. 8. Thus, this 
solution has not been programmed in the VAST code; however, it will be considered as an 
alternate procedure to compare various pressure solutions. Its deficiencies for the 
damping seaij analysis are the inefficient iterative proceaure for time dependent, three- 
dimensional solutions and the requirement of Neuman boundary conditions. 

The artificial compressibility method, Eq. (6) was programmed with the VAST code for 
evaluation. Results for low Reynolds number flows presented in Ref. 6 were good. Since 
an Identical equation was used in the VAST code, only high Reynolds number cases were 
examined, the majority of time with liquid (water) as the working fluid. For pipe flow, 
the velocities generally developed well, but not completely. However, pressures were 
never correct within the pipe and, in fact, showed more of a lateral variation than the 
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requisite axial dependence. Results for B = 1650 are shown In Fig, 2, This value 
showed the best results of six different values between 1.0 and 8100. Because of the 
dynamic differencing procedure of the VAST code, limits for 8 , which cccording to Ref. 
5 are dependent on the numerical scheme, could not be accurately developed. Hence a 
trial and error procedure was necessary. The parameter 8 also appeared in the varia- 
tional functional. This led to a low level of confidence in the ability of the functional to 
accurately dictate the differencing schemes. Due to these problems, it was deemed that 
the artificial compressibilities pressure solution should not be used If another procedure 
was found to be adequate. 

To attempt to overcome the problems of existing pressure solution methods, the 
procedure represented by Eq. (II) was developed. Results for both inviscid and viscous 
flow have been obtained. For a cylinder in a crossf'ow (half-plane) the Inviscid solution 
is shown in Fig. 3. This figure is for a comparatively fine grid of 1,008 elements; a grid 
consisting of 224 elements was also run. The results for the fine grid were closer to the 
analytical solution, although both cases were qualitatively good. The fine grid was run in 
an effort to produce better results just downstream of the cylinder. Although the 
streamlines can be shown to be good, the magnitudes of velocities and the pressures are 
seen to be less than expected. This is beOved to be partially due to the initial 
conditions, which were the assumption of the uniform freestream velocity and pressure 
everywhere, including the node immediately behind the cylinder. 

Several sets of boundary conditions were also investigated for this problem. The 
conditions used to produce the results of Fig. 3 were: slip walls for the cylinder, 
centerline, and top of the grid, static pressure on the downstream boundary, and 
velocities (magnitude and direction) at the upstream boundary. One boundary condition 
which has not yet been investigated, but which could possibly solve the downstream flow 
discrepancy, is the use of constant pressure gradient normal to the cylinder surface. This 
condition, however, would eliminate one advantage of this procedure over a Poisson 
equation. 

Viscous flow over the cylinder was also studied with the course grid. Reynolds numbers 
of i.6 x 10^ and 1.6 x 10^ were run with laminar viscosity values. Results of both cases 
were good, with those of the higher R^ shown in Fig. 4. The separation and recirculation 
behind the cylinder are quantitatively very good. Pressure distribution Is also correct. 
Boundary conditions similar to those of the inviscid case were used; however, surface 
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Fig. 2a Pipe Flow Calculated Using Artificial Compressibility 
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force which was suggested in Ref. 12 was not used. It Is expected that this condition 
would improve the already good solution. 

Pipe flow problems, both two-dimensional and axlsymmetrlc, have also been studied. 
Results for air flow (axlsymmetrlc) at Rfvj = 1000 are shown in Fig. 5. The velocity 
profile developed well, though not completely, and the pressure is seen to be essentially a 
function of axial distance alone, as it should be. The predicted pressure drop of 
approximately 24 psf compares well to the analytical value of 21.4 psf for this problem. 
Tills problem has not reached steady state, and further integration should improve the 
accuracy. 

With water as the fluid, pipe flow results were not as encouraging. For a Reynolds 
number of 2.78 x 10^, the velocities, shown In Fig. 6, are developing well, though 
slowly. The pressure solution Is poor, however. The distribution varies axially, as 
opposed to the artificial compressibility, but pressure decreases In the upstream 
direction and even takes on forge negative values. Although this problem is not fully 
understood, the P r c^ term of Eq. (II) completely dominates for this flow (10^ times as 
large as the P and P r terms). The large velocity derivatives produced at the early time 
steps then effectively swamp the pressure solution. Methods of alleviating this problem 
other than developing realistic initial profiles with small gradients are under study. 
Boundary conditions are also being investigated as a potential source of error. The 
surface force equation (slightly modified) of Ref. (12) was used as the downstream 
boundary condition. Mass flow and normal pressure gradient boundary conditions could 
also be required for the pipe problem. However, the slowly developing solution may 
simply not yet be close enough to steady-state; this possibility is also being investigated 
further. 

In summary, the pressure solution given by Eq. (II) shows promise, but it has not been 
developed sufficiently. However, it does appear to be the best method if an accurate 
pressure field is desired. The only alternative is a Poisson equation. Long run times 
seem to be associated with both methods, as the small CFL produced by the high acoustic 
speeds of liquids requires many integration steps for the explicit pressure solution. One 
important result of the success achieved with Eq, (II) for the pressure was the 
performance of the entropy functional as a variational principle. Although the principle 
was hypothesized and used for compressible, ideal gases, no stability problems were 
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encountered for the Incompressible cases run. This Is certainly a factor In pursulnq Eq. 
(I I) as the best pressure solution. 


CLOSURE 

Once the pressure equation has been sufficiently developed, the main obstacle to the 
damping seals analysis will have been cleared. Eccentric operation and three-dimenslnal 
solutions are straightforward extensions of the VAST code. Roughness effects can be 
handled empirically by the turbulence model developed as Eqs. (19) and (20). Transport 
equation models are also under evaluation. Damping seals narametrlc cases can then be 
studied. 
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